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ABSTRACT 

We excite an epicyclic motion, whose amplitude depends on the vertical position, z, 
in a simulation of a turbulent accretion disc. An epicyclic motion of this kind may be 
caused by a warping of the disc. By studying how the epicyclic motion decays we can 
obtain information about the interaction between the warp and the disc turbulence. 
A high amplitude epicyclic motion decays first by exciting inertial waves through a 
parametric instability, but its subsequent exponential damping may be reproduced by 
a turbulent viscosity. We estimate the effective viscosity parameter, ctv, pertaining 
to such a vertical shear. We also gain new information on the properties of the disc 
turbulence in general, and measure the usual viscosity parameter, ah, pertaining to a 
horizontal (Keplerian) shear. We find that, as is often assumed in theoretical studies, 
Ofv is approximately equal to ah and both are much less than unity, for the field 
strengths achieved in our local box calculations of turbulence. In view of the smallness 
(~ 0.01) of av and ah we conclude that for (3 = Pgas/Pmag ^ 10 the timescale for 
diffusion or damping of a warp is much shorter than the usual viscous timescale. 
Finally, we review the astrophysical implications. 
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1 INTRODUCTION 

Warped accretion discs appear in many astrophysical sys- 
tems. A well known case is the X-ray binary Her X-1, in 
which a precessing warped disc is understood to be periodi- 
cally covering our line of sight to the neutron star, resulting 
in a 35-day periodicity in the X-ray emission (Tananbaum et 
al. 1972; Katz 1973; Roberts 1974). A similar phenomenon 
is believed to occur in a number of other X-ray binaries. In 
recent years the active galaxy NGC 4258 has received much 
attention as a warp in the accretion disc has been made 
visible by a maser source (Miyoshi et al. 1995). 

A warp may appear in an accretion disc in response to 
an external perturber such as a binary companion, but it is 
also poss ible t hat the disc may produce a warp on its own. 
Pringle (1996) showed that the radiation pressure from the 



central radiation source may produce a warp in the outer 
disc. In a related mechanism the irradiation can drive an 
outflow from the disc. The force of the wind may then in 
a similar way excite a warp in the disc (Schandl & Meyer 
1994). 

Locally, one of the effects of a warp is to induce an 
epicyclic motion whose amplitude varies linearly with dis- 
tance from the midplane of the disc. This motion is driven 
near resonance in a Keplerian disc, and its amplitude and 
phase are critical in determining the evolution of the warp 
(Papaloizou & Pringle 1983; Papaloizou & Lin 1995). De- 
pending on the strength of the dissipative process the warp 
may either behave as a propagating bending wave or evolve 
diffusively. In the latter case the amplitude of the epicyclic 
motion is determined by the dissipative process. 

A (possibly) related dissipative process is responsible 
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for driving the inflow and heating the disc by transporting 
angular momentum outwards. From a theoretical point of 
vie w this transport has been described in terms of a viscos- 
ity ( shakura & Sunyaev 1973 ) , but the source of the viscosity 
remained uncertain for a long time. It was clear from the be- 
ginning that molecular viscosity would be insufficient, so one 
appealed to some form of anomalous viscosity presumably 
produced by turbulence in the accretion disc. However the 
cause of the turbulence could not be found as the Keplerian 
rotation is hydrodynamically stable according to Rayleigh's 
criterion. 

Eventually Balbus 



Hawley (1991) discovered that 



the Keplerian flow becomes unstable in the presence of a 
magnetic field. This magnetic sh earing instability had al- 
ready been described by Velikhov ( 1959 ) and Chandrasekhar 
(1960), but it had not been thought applicable in the con- 
text of accretion discs before. Several numerical simulations 
(e.g. Hawley, Gammie & Balbus 1995, Matsumoto & Tajima 
1995, Brandenburg et al. 1995, Stone et al. 1996) have 
demonstrated how this instability generates turbulence in 
a Keplerian shear flow. 

The most important result of these simulations has been 
to demonstrate that the Maxwell and Reynolds stresses that 
the turbulence generates will transport angular momentum 
outwards, thus driving the accretion. The energy source of 
the turbulence is the Keplerian shear flow, from which the 
magnetic field taps energy. This energy is then partially dis- 
sipated due to Ohmic diffusion, but an equal amount of en- 
ergy is spent on exciting the turbulent motions. The turbu- 
lent motions on the other hand give rise to a dynamo, which 
sustains the magnetic field required by the Balbus-Hawley 
instability (Brandenburg et al. 1995, Hawley, Gammie & 
Balbus 1996). 

In general the energy of the magnetic field is an order of 
magnitude larger than the energy of the turbulent velocities, 
but almost all of the magnetic energy is associated with 
the toroidal magnetic field, and the poloidal magnetic field 
components are comparable to the turbulent velocities. 

So far none of the simulations has addressed the ques- 
tion of how the turbulence responds to external perturba- 
tions or systematic motions that are more complex than a 
Keplerian shear flow. The purpose of this paper is to begin 
such an investigation by studying how the turbulence inter- 
acts with an imposed shearing epicyclic motion of the type 
found in a warped disc. 

We start this paper by describing the shearing-box 
approximation of magnetohydrodynamics and summarizing 
the properties of the epicyclic motion in a shearing box in 
Sect. 2. Section 3 is then a description of our simulations of 
an epicyclic motion. The results of the simulations are then 
described in Sect. 4 and briefly summarized in Sect. 5. 



where p is the pressure, p the density, and gz = —GMz/Rq 
the vertical component of the gravity with G the gravita- 
tional constant, M the mass of the accreting star, and Ro 
the radial distance from the star. For simplicity we assume 
that the disc material is initially isothermal, and is a perfect 
gas, so that p = c^p, where Cs is the isothermal sound speed, 
which is initially constant. The density distribution is then 



P = poe 



where the Gaussian scale height, H , is given by 
2c! Rl 



H 



2.2 



GM ' 

Epicyclic motion in the shearing box 
approximation 



(2) 



(3) 



In the shearing box approximation a small part of the accre- 
tion disc is represented by a Cartesian box which is rotating 
at the Keplerian angular velocity Qo — ^ GM/Rq. The box 
uses the coordinates {x, y, z) for the radial, azimuthal and 
vertical directions, respectively. The Keplerian shear flow 
within the box is u^^ = —^flox, and we solve for the devi- 
ations from the shear flow exclusively. The magnetohydro- 
dynamic (MHD) equations may then be written 



Vp 

vt " 

Vu 

VB 

Vt 



-V ■ (pu) 



(4) 



- (u ■ V) u+g+f (u)-ivp+ij X B+iv- (2upS) , (5) 



: Vx (u X B) 



Vxj7/i()J, 



(6) 



^ ^ _ (u . V) e-^V ■ u+-\/ixp\7e)+2t.S'+'-l^J^+Q, (7) 
Vt p p p 

where T>/T>t = d/dt + u^^d/dy includes the advection by 
the shear flow, p is the density, u the deviation from the Ke- 
plerian shear flow, p the pressure, f(u) — Qo{2uy,—^Ux,0) 
the inertial force, B the magnetic field, J = V x B//io the 
current, po the permeability of free space, f the viscosity, 
Sij ~ + ~ ^SijUk^k) the trace-free rate of strain 

tensor, 77 the magnetic diffusivity, e the internal energy, x the 
thermal conductivity, and Q is a cooling function. The ra- 
dial component of the gravity cancels against the centrifugal 
force, and the remaining vertical component is g = —Qfizz. 
We adopt the equation of state for an ideal gas, p — (7— l)/9e. 

When the horizontal components of the momentum 
equation ^ are averaged over horizontal layers (an oper- 
ation denoted by angle brackets), we obtain 



^/ \ on/ \ ^/ \^^ / 

-^ApUx) = 2\lo{pUy) - TT-ipUxUz) + { 

Ot OZ OZ \ po 



(8) 



2 MATHEMATICAL FORMULATION 

2.1 The local structure of a steady disc 

For the intentions of this paper it is sufficient to use a simple 
model of the vertical structure of a geometrically thin accre- 
tion disc. The disc is initially in hydrostatic equilibrium. 



dp 

OZ 



(1) 



-^(pUy) = --ilo{pUx) - ^{pUyUz) + — ( 

Ot 2 OZ OZ \ po 



(9) 



The explicit viscosity, which is very small, has been ne- 
glected here. These equations contain vertical derivatives of 
components of the turbulent Reynolds and Maxwell stress 
tensors, distinct from the 2; y- components that drive the ac- 
cretion. 

We initially neglect the turbulent stresses and obtain 
the solution 
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{pux) = pq{z)uo (z) cos{flot), 

1 



(pUy) 



-po(z)uo (z) sin(r2oi), 



(10) 
(11) 



which describes an epicychc motion. Here po{z) is the ini- 
tial density profile. The initial velocity amplitude uo is an 
arbitrary function of z. For the simulations in this pa- 
per we will take uo{z) oc sm{kz), where k — n/L^, and 
— \Lz < z < ^Lz is the vertical extent of our shearing 
box. This velocity profile is compatible with the stress-free 
boundary conditions that we employ in our numerical sim- 
ulations, and gives a fair representation of a linear profile 
close to the midplane of the disc. 

The kinetic energy of the epicyclic motion is not con- 
served, but the square of the epicyclic momentum 



E{z,t) = - {pu^f +2{puyf , 



(12) 



is conserved in the absence of turbulent stresses. By multi- 
plying Eq. ^ by {pu^}, and Eq. (^ by 4{puy) we obtain 

dE_ 
dt 

where 



Fu + Fb, 



(13) 



d d 

Fu = ~ {py-x) {pUxUz) — 4 {pUy) — {pUyUz) 

and 



Fb = (pUj:) 



d I B^Bz 
dz 



po 



dz 



ByBz 

po 



(14) 



(15) 



represent the 'rates of working' of the Reynolds and Maxwell 
stresses, respectively, on the epicyclic oscillator. We may ex- 
pect that both Fu and Fb are negative, but by measuring 
them in the simulation we may determine the relative impor- 
tance of the Reynolds and Maxwell stresses in damping the 
epicyclic motion. We will also refer to an epicyclic velocity 
amplitude 



u = y/{u^y + Huyy. 



2.3 Theoretical expectations 



(16) 



The detailed fluid dynamics of a warped accreti on dis c has 
been discussed b y, e.g ., Papaloizou fc Pri ngle ( 1983 ), Pa- 
paloizou & Lin ( 1995 ) and Ogilvie (199£). The dominant 
motion is circular Keplerian motion, but the orbital plane 
varies continuously with radius r and time t. This may con- 
veniently be described by the tilt vector £{r,t), which is a 
unit ve ctor parallel to the local angular momentum of the 




disc aniulus at radius r. A dimensionless measure of the 



Figure 1. Owing to the stratification there appear horizontal 
pressure gradients in the warped disc. These gradients excite the 
epicyclic motion (arrows indicate forces, not velocities) 



not yet possible to study simultaneously both the small- 
scale turbulence and the global dynamics of the accretion 
disc in a numerical simulation. One of the goals of this pa- 
per is to test the validity of this hypothesis by comparing 
the predictions of the viscous model, as summarized be- 
low, with the results of the numerical model. We generalize 
the viscosity prescription by allowing, in a simple way, for 
the possibility that the effective viscosity is anisotropic (cf. 
Terquem 1998). The parameter ah pertaining to 'horizontal' 
shear (i.e. horizontal-horizontal components of the rate-of- 
strain tensor, such as the Keplerian shear) may be different 
from the parameter pertaining to 'vertical' shear (i.e. 
horizontal-vertical components of the rate-of-strain tensor, 
such as the shearing epicyclic motion). 

Owing to the pressure stratification, resulting from the 
vertical hydrostatic equilibrium, in a warped disc there are 
strong horizontal pressure gradients (Fig. ^) , which generate 
horizontal accelerations of order AQqz, that oscillate at the 
local orbital frequency, as viewed in a frame co- rotating with 
the fluid. In a Keplerian disc the frequency of the horizontal 
pressure gradients coincides with the natural frequency of 
the resulting epicyclic motion, and a resonance occurs. The 
amplitude of the resulting epicyclic motion depends on the 
amount of dissipation present. At low viscosities, Qv ^ H/r, 
the amplitude is limited by the coupling of the epicyclic mo- 
tion to the vertical motion in a propagating bending wave, 
which transports energy away. At higher viscosities the am- 
plitude is limited by a balance between the forcing and the 
viscous dissipation, and the warp evolves diffusively ( |^a-| 
paloizou & Pringle 198o). When H/r ay <S 1 the ampli- 



amplitude of the warp is then A — \d£/d\nr\. 

In the absence of a detailed understanding of the tur- 
bulent stresses in an accretion disc, it is invariably assumed 
that the turbulence acts as an isotropic effective viscosity in 
the sense of the Navier-Stokes equation. In such an approach 
the dynamic viscosity is often parametrized as 



tude of the epicyclic motion is 



(18) 



fi = ap/Qo, 



(17) 



where a is a dimensionless parameter (Shakura & Sunyaev 
1973). Although it is now possible to simulate the local tur- 
bulence in an accretion disc, this form of phenomenological 
description of the turbulent stress is still valuable as it is 



The resulting hydrodynamic stresses pUxUz and pUyUz (over- 
bars denote averages over the orbital timescale) , which tend 
to flatten out the disc, are also proportional to a^^ A and 
therefore dominate over the stresses oc a^A due to small- 
scale turbulent motions, which would have the same effect. 
It is for this reason that the timescale for flattening a warped 
disc is anomalously short compared to the usual viscous 
timescale, by a factor of approximately 2ahav. (For more 
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details, see Papaloizou & Pringle 1983 and Ogilvie 1999. In 
these papers it was assumed that ah = Ov) 

We aim to measure both Qh and Ov in the simulations. 
The first may be obtained through the relation 



BxBy \ 3 , , 

pu^Uy ) = 7Tah<p)^ 

MO /v 2 



(19) 



which follows by identifying the total turbulent xy-stiess 
with the effective viscous xy-stress resulting from the vis- 
cosity acting on the Keplerian shear (note that the def- 
inition of Oh differs from that of ass in Brandenburg et al. 
1995 by a factor \/2). Here the average is over the entire 
computational volume. 

The coefficient av may be obtained by measuring the 
damping time of the epicyclic motion. In the shearing-box 
approximation the horizontal components of the Navier- 
Stokes equation for a free epicyclic motion decaying under 
the action of viscosity are 



dux „^ 19/ OvP dux \ 

at p oz \ uo oz ) 

duy 1^ 19 fa.vpduy\ 

— VlnU-r H — 

dt 2 " " p 92 V no dz J 



(20) 
(21) 



Under the assumptions that av and dzU are independent of 
z and that the disc is vertically in hydrostatic equilibrium, 
these equations have the exact solution 



Ciloze 
1 



cos(r2oi), 
Ciloz e"*'''" sin(f2of). 



where C is a dimensionless constant and 
1 



avQo 



(22) 
(23) 

(24) 



is the damping time. Admitt edly it is already believed that 
Qh is not independent of z (Brandenburg et al. 1996) and 
therefore Ov may not be either. Also our velocity profile 
is not exactly proportional to z. However, it is in fact the 
damping time that matters for the application to a warped 
disc, and the solution that we describe above is in a sense the 
fundamental mode of the epicyclic shear flow in a warped 
disc. 

It might be argued that the decaying epicyclic motion 
is fundamentally impossible in the presence of MHD tur- 
bulence. After all, the magnetorotational instability works 
because a magnetic coupling between two fluid elements ex- 
ecuting an epicyclic motion allows an exchange of angular 
momentum that destabilizes the motion. However, we argue 
that the epicyclic motion must be re-stabilized in the nonlin- 
ear turbulent state. Otherwise epicyclic motions, which are 
continuously and randomly forced by the turbulence, would 
grow indefinitely ( Balbus fc Hawley 199^ ). In fact they last 
only a few orbits, as discussed in Section 3.2 below. More- 
over, our numerical results demonstrate without any doubt 
that the shearing epicyclic motion is possible, and does de- 
cay in the presence of MHD turbulence. 

A further theoretical expectation is as follows. In an 
inviscid disc, the epicyclic motion can decay by exciting 
inertial waves through a parametric instability (Gammie, 
Goodman & Ogilvie 2000, submitted). In the optimal case, 
the signature of these waves is motion at 30° to the vertical. 



Table 1. Specification of numerical simulations: the number of 
grid points are given by X Ny X , and the amplitude of the 
initial velocity perturbation is uq. Note that all Runs except Run 
4 starts from a snapshot of a previous simulation of turbulence 
in an accretion disc. Run 4 starts from a state with no turbulent 
motions 



Run 


Nj; X Ny X Nz 


Uo 





31 X 63 X 63 


0.0 


1 


31 X 63 X 63 


0.011 


lb 


63 X 127 X 127 


0.011 


2 


31 X 63 X 63 


0.095 


3 


63 X 127 X 127 


0.095 


4 


63 X 127 X 127 


0.095 



while the wave vector is inclined at 60° to the vertical, 
characteristic local growth rate of the instability is 



The 



7 = 



3\/3 
16 



9l 



dz 



(25) 



This instability can lead to a rapid damping of a warp, 
but may be somewhat delicate as it relies on properties 
of the inertial- wave spectrum. It is important to determine 
whether it occurs in the presence of MHD turbulence. 



3 NUMERICAL SIMULATIONS 
3.1 Computational method 



We use the code by Nordlund & Stein (|1990|) with the 
modif ications that were described by Brandenburg et al. 
( [1995| ). The code solves the MHD equations for In p, u, e and 
the vector potential A, which gives the magnetic fleld via 
B = V X A. For the (radial) azimuthal boundaries we use 
(sliding-) periodic boundary conditions. The vertical bound- 
aries are assumed to be impenetrable and stress- free. Unlike 
our earlier studies, we now adopt perfectly conducting ver- 
tical boundary conditions for the magnetic fleld. Thus we 
have 



9m, 



^ = = 0, 
oz 



dBy 



B, = 0. 



(26) 



(27) 



dux 
and 

dBx ^ 

dz dz 

We choose units such that H = GM — 1. Density is 
normalized so that initially p = 1 at the midplane, and we 
measure the magnetic fleld strength in velocity units, which 
allows us to set po — 1. The disc may be considered to 
be thin by the assumptions of our model, and the results 
will thus not depend on the value of _Ro. We choose to set 
Ro = 10 in our units, which gives the orbital period Tq = 
2tv/Q,o = 199, and the mean internal energy e = 7.4 10~ . 

= 1 : 27r : 4, where x and 
respectively, and y goes 



The size of the box is : Ly : 
z vary between ±^Lx and ±^Lz 



To stop the box from heating up during the simulation we 
introduce the cooling function 



= — o-cooi - 



eo) , 



(28) 



where CTcooI is the cooling rate, which typically corresponds 
to a timescale of 1.5 orbital periods, and eo is the internal 
energy of an isothermal disc. 



© 0000 RAS, MNRAS 000, 000-000 



Vertical shear flow and turbulence 5 



Although the Balbus-Hawley instability appears readily 
in a numerical simulation, experience has taught us that the 
initial conditions must be chosen carefully in order for the 
instability to develop into sustained turbulence, unless the 
initial magnetic field has a net flux. In particular, the initial 
field strength should be chosen such that the Alfven speed 
is close to the sound speed. Otherwise the field becomes too 
weak before the dynamo effect sets in. Experience has shown 
that after about 50 orbital periods the simulations are in- 
dependent of the details of the initial conditions, which is 
typical of turbulence in general. To save time we started 
the simulat ions in this paper fr om a snapshot of a previous 
simulation (Brandenburg 1999). The origin of this sn apsho t 
goes back to the simulations by Brandenburg et al. ( 1995 ). 
In those simulations we started from a magnetic field of the 
form zBo sin{2nx/ Lx) (at that time we employed boundary 
conditions that constrained the magnetic field to be vertical 
on the upper and lower boundaries). Bo was chosen such 
that (3 = 2 flop /B^ = 100 on average in the shearing box. 
A snapshot from an evolved stage of one of these simula- 
tions was later relaxed for about 55 orbital periods to fit the 
perfectly conducting upper and lower boundaries of Bran- 
denburg (1999). For the purposes of this paper a snapshot 
from the new simulation was modified in the following way. 
For every horizontal layer in the snapshot we subtract the 
mean horizontal velocity and then add a net radial flow of 
the form 



Ux ~ uo sin 



(29) 



Two things that should be kept in mind here are that the 
previous simulations have been run long enough that the 
current snapshot has lost its memory of its original initial 
conditions, and that in spite of all the modiflcations there is 
still no net magnetic flux in the shearing box. The number of 
grid points and uq for the different runs are given in Table |l| 
Mo should be compared to the adiabatic sound speed which 
is 0.029. We include one run, Run 0, in which we do not 
excite an epicyclic motion, as a reference. 



3.2 Results 

We start by looking at Run 0. In this case we have not mod- 
ifled the velocity fleld, that is Run represents the typical 
properties of the turbulence. The turbulence is fully devel- 
oped at the beginning of our study meaning that the shear- 
ing box does not remember its initial state any longer. We 
plot the vertical distributions of the magnetic and turbulent 
kinetic energies in Fig. |[ The kinetic energy is indepen- 
dent of the vertical coordinate z while the magnetic energy 
has a minimum close to the midplane. It is a typical prop- 
erty of all our simulations that the magnetic energy is up 
to an order of magnitude larger than the turbulent kinetic 
energy, but the magnetic energy is still an order of magni- 
tude smaller than the internal energy, whose mean density 
as stated above is 7.410^*. A property of both the mag- 
netic and turbulent kinetic energies is that they are com- 
pletely dominated by the contributions from the azimuthal 
components of the magnetic field and turbulent velocity, re- 
spectively. The other components of the turbulent velocity 
and magnetic field make about equal contributions to the 
energy density far from the midplane of the disc. 
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Figure 2. Top: The magnetic (solid line) and kinetic (dashed 
line) turbulent energies as a function of z at the beginning of 
Run 0. The energies are completely dominated by the contribu- 
tions from the ^/-components of the magnetic field, and velocity, 
respectively. Bottom: (dashed line), (dot-dashed line), 

(ni^l2 (triple-dotted dashed line) and pu\l2 (dotted line) as a 
function of z 
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Figure 3. Top: {ux) (solid line) and 2(iiy) (dashed line) as a 
function of time at 2; = 1.55 for Run 0. The plot shows that 
epicyclic motions of amplitude ~ 0.005 lasting for a couple of 
orbits appear spontaneously in the disc. Bottom: m as a function 
of z at t = 66.7 orbital periods 
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Figure 4. {ux) (solid line) and 2{uy) (dashed line) as a function 
of time at z = 1.17 (top) and z = —0.25 (bottom) of Run lb. The 
damping timescale in the lower plot during the interval 63 to 68 
orbital periods is 15.5 orbital periods 
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Figure 5. (E)y, the square of the epicyclic momentum, as a 
function of time for Run lb. The dashed line is a fit with an 
exponential to (E)v- The e-folding timescale of the exponential 
is 17.8 orbital periods 

The turbulence generates epicyclic motions with ampli- 
tudes ~ 0.004 lasting for a couple of orbital periods (Fig. 

The epicyclic velocity amplitude u (see Eq. p^ ) is peaked 
towards the surfaces of the accretion disc (cf. Fig. ^, bot- 
tom). Overall these motions complicate the analysis of the 
rest of our numerical simulations, and force us to excite mo- 
tions with amplitudes significantly larger than that of the 
motions produced by the turbulence. 

Figure ^ shows the mean horizontal motion of Run lb. 
At the start of the simulation we added a radial velocity 
of amplitude 0.011. Far from the midplane the imposed 



^2 ^1 1 2 

z 

Figure 6. il as a function of z at t = 55.8 To (solid line), 57.1 To 
(dashed Une), 58.4 To (dotted line), and 60.9 To (dot-dashed line) 
for Run 3 

epicyclic motion is comparable to that generated by the tur- 
bulence, and it becomes virtually impossible to identify a 
phase of exponential decay. Closer to the midplane the im- 
posed epicyclic motion is initially unaffected by the turbu- 
lence. After seven or eight orbital periods a damping sets in, 
but the net damping before the turbulence starts to reinforce 
the epicyclic motion is less than a factor of two. By looking 
at the volume average of the square of the epicyclic mo- 
mentum, {E)y, (Fig. n) we find a clear trend. The epicyclic 
motion is obviously damped, and by fitting an exponential 
we obtain an e-folding timescale of 17.8 orbital periods for 
{E)v, corresponding to a timescale of 35.6 orbital periods 
for the momentum, but the damping is not described ac- 
curately by an exponential. Initially the damping is much 
slower than the fitted exponential, while at the end it is 
faster. 

The analysis becomes more straightforward in Run 3, 
where we excite a motion of amplitude 0.095. There is then a 
sufficient dynamic range between the epicyclic and turbulent 
motions. We show the vertical variation of u at four different 
times in Fig. and plot it as a function of time at three 
different heights; z = 1.52, z = 0.89 and z = 0.44, in Fig. 0. 
Figure ^ shows that the damping sets in first at the surfaces, 
while for \z\ < 1 there is essentially no damping during 
the first two orbital periods (Fig. There is then a brief 
period of rapid damping between t = 58 To and t = 60 Tq 
throughout the box, especially for small \z\ where u may 
drop by a factor 2. This is followed by a period of exponential 
decay, but after t = 75To it becomes difficult to follow the 
epicyclic motion, as the influence of the random turbulence 
on u becomes significant, in particular close to the midplane, 
where u is anyway small. We estimate the damping time, 
r — (dlnu/dt)~^ by fitting exponentials to u in the interval 
60To <t< 75To. Averaged over the box we get r = 25±8ro, 
which corresponds to Qv = 0.006 ± 0.002 according to Eq. 

(§)• 

We may determine the influence of the Maxwell and 
Reynolds stresses on the shear flow by plotting Fu and Fb 
as functions of time (Fig. ^ . The sharp peak in F„ coincides 
with the parametric decay of the epicyclic motion to inertial 
waves (see Sect. ^^). As this is a hydrodynamic process it is 
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Figure 7. The amplitude of the epicychc motion in Run 3, u, 
as a function of t on three horizontal planes: z = 1.52 (solid 
line), 2 = 0.89 (dashed line) and z = 0.44 (dot-dashed line). 
The straight lines are exponential functions that have been fitted 
for the interval 60To < t < 75To. The e-folding timescales of 
the exponentials starting from the top are 19.4ro, 22.iro and 
20.2To, respectively. The epicyclic motion was added to the box 
at t = 55.7ro 
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Figure 8. —{Fu)v (solid line) and — (Ffl)v (dashed line) as a 
function of time for Run 3 

not surprising that Fu > Fb, but it was not expected that 
Fu would remain the dominating efltect at later times, since 
at the same time it is the Maxwell stress that is driving the 
radial accretion flow. 

The accretion itself is driven by the (pu^Uy — BxBy)- 
stresses. We plot the moving time averages of the vertical 
average of the accretion-driving stress of Run 3 in Fig. ^ 
and as a comparison that of Run in Fig. In the be- 
ginning of Run 3 the Reynolds stress is modulated on a 
timescale of half the orbital period. This modulation is an 
artifact of the damping of the epicyclic motion and dies out 
with time. The Maxwell stress becomes significantly stronger 
than the Reynolds stress once the epicyclic motion has van- 
ished. Later the stresses vary in phase with each other, like 
they do all the time in Run 0. The main difference between 
Run and the end of Run 3 is that the stresses are 2-3 times 



Figure 9. —{BxBy)\r (solid line) and (pUxUy)^ of Run 3 
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Figure 10. ~{BxBy)Y (solid line) and {pUxUy)\r of Run 
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Figure 11. a^^ is calculated by dividing —{BxBy)^ (solid line) 
and (pUxUy)Y (dashed line) of Run 3 with ^(p)v 



© 0000 RAS, MNRAS 000, 000-000 



8 U. Torkelsson et al. 



0.010 



0.005 



m 0.000 



-0.005 - 



-0.010 




0.10 



60 65 70 75 80 85 
t (orbital periods) 



23 



0.01 



1.52 



0.89 



0.44 



2 4 6 

t (orbital periods) 



Figure 12. By averaged over the upper half (solid line) and lower 
half (dotted line) of the box of Run 3 



larger in Run 3. We calculate ah for Run 3 by dividing the 
stresses by |{p)v (Fig- 0)- One should take note that the 
simulations in this paper are too short to derive Qh with high 
statistical significance (cf. Brandenburg et al. 1995) , but our 
results do show that Oh varies in phase with the stress. In 
other words the pressure variations are smaller (the pres- 
sure increases by 50% during the course of the simulation) 
than the stress variations, which is not what we expect from 
Eq. dig]). The lack of a correlation between the stress and 
the pressure is even more evident in Run 0, in which the 

pressure never varies with m ore than 5%. 

In our previous work (Brandenburg et al. 1995) we 



found that the toroidal magnetic flux reversed its direction 
about every 30 orbital periods. With the perfect conductor 
boundary conditions that we have assumed in this paper 
such field reversals are not allowed, since the boundary con- 
ditions conserve the toroidal magnetic flux (e.g. Branden- 
burg 1999). In all of our simulations we made sure that the 
toroidal magnetic flux was 0. However the azimuthal mag- 
netic fleld organized itself in such a way that it is largely 
antisymmetric with respect to the midplane. That is, con- 
sidering separately the upper and lower halves of the box, 
we may still find significant toroidal fluxes, but with op- 
posite directions. These fluxes may still be reversed by the 
turbulent dynamo as we found in Run 3 (Fig. 



3.3 Dependence on uq and on the resolution 

There is some evidence from comparing Run lb and Run 3 
that the damping timescale decreases as mq is increased, but 
it is more difficult to study the damping of the epicyclic mo- 
tion for a smaller uq, as the turbulence may excite epicyclic 
motions on its own. The randomly excited motions swamp 
the epicyclic motion that we are studying. The quantitative 
results from Runs 1 and lb are therefore more uncertain. In 
other respects the turbulent stresses of Runs 1 and lb are 
more similar to those of Run than to those of Run 3. 

On the other hand there are significant differences be- 
tween Runs 2 and 3, which differ only in terms of the grid 
resolution. The magnetic field decays rapidly in Run 2, and 
only the toroidal fleld recovers towards the end of the sim- 



Figure 13. The amplitude of the epicyclic motion of Run 4, u, 
as a function of t on three horizontal planes: z = 1.52 (solid 
line), 2 = 0.89 (dashed line), and z = 0.44 (dot-dashed line). 
The straight lines are exponential functions that have been fitted 
for the interval < t < 5.5To. The e-folding timescales of the 
exponentials starting from the top are 640To , 4 700To and 2 SOOTq , 
respectively 



ulation. In the absence of a poloidal magnetic fleld there is 
no magnetic stress, and therefore the disc cannot extract 
energy from the shear flow. Consequently there is no tur- 
bulent heating in Run 2, and the disc settles down to an 
isothermal state. Apparently the turbulence is killed by the 
numerical diffusion in Run 2. This demonstrates that the 
minimal resolution which is required in the simulation de- 
pends on the amplitude uq. A simulation with an imposed 
velocity uq with an amplitude signiflcantly larger than that 
of the turbulence requires a flner resolution than a simu- 
lation of undisturbed turbulence. Run 2 failed because the 
imposed velocity fleld in combination with the limited res- 
olution generated a numerical diffusion large enough to kill 
the turbulence. This was not a problem in Run 1, where the 
amplitude uq is much smaller. 

To check the influence of numerical diffusion on Run 
3 we imposed the same velocity proflle as in Run 3 on a 
shearing box without any turbulence (Run 4). We plot the 
evolution of u in Fig. y_3|. We fltted exponential functions 
to u for the interval to 5.5 orbital periods. The shortest 
e-folding time we obtained over this interval was 640 orbital 
periods for z = 1.52. Closer to the midplane the e-folding 
timescale was several thousand orbital periods. These num- 
bers are indicative of the effect of numerical diffusion in the 
simulations. The strong damping at t = 6 orbital periods 
is not due to numerical diffusion, rather it i s caused by a 
physical parametric instability (see Sect. 4.1) that, in this 
case, has been triggered by numerical noise. 



4 DISCUSSION 

For the purposes of the following discussion we now summa- 
rize our main results: 

• In Runs 3 and 4 the initial damping of the epicyclic 
motion is cau sed by its parametric decay to inertial waves 
(see Sect. O). 
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Figure 14. (u'^)v (top) and {u'^)v (bottom) for a two- 
dimensional simulation of parametric decay. The decay of the 
epicyclic motion excites a vertical motion 
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Figure 15. The density (p) (top) and internal energy (e) (bot- 
tom) as functions of z at the start of the two-dimensional simu- 
lation (solid line) and after 8.97 orbital periods (dashed line) 



• Apart from this, the epicyclic motion experiences ap- 
proximately exponential damping through interaction with 
the turbulence. The e-folding timescale is about 25 orbital 
periods, which may be interpreted as — 0.006. 

• Qh, which describes the accretion-driving stress 



{pUx 



BxBy) is of comparable size. 



4.1 Parametric decay to inertial waves 



Gammie et al. (200C ) predicted the occurrence of a paramet- 
ric instability in epicyclic shear flows. The shear flow should 
excite pairs of inertial waves that propagate at roughly a 30° 
angle to the vertical, and involve vertical as well as horizon- 
tal motions. To elucidate the dynamics of the parametric 
instability we make use of a two-dimensional hydrodynamic 
simulation of an epicyclic shear flow. Our two-dimensional 
xz-plane has the same extension as in the previous three- 
dimensional simulations, but we are now using 128 x 255 
grid points. The initial state is a stratified Keplerian ac- 
cretion disc to which we have added a radial motion with 
amplitude 0.095, and a small random perturbation of the 
pressure. Initially the amplitude of the epicyclic motion, u, 
is constant, but the damping sets in suddenly after 5 orbital 
periods (Fig. |lj). Over four orbital periods the amplitude 
of the epicyclic motion decreases by 30%, after which the 
damping becomes weaker again. This is clearly not an expo- 
nential damping. At the same time as the epicyclic motion 
is damped the vertical velocity starts to grow. The iner- 
tial waves themselves are not capable of transporting mass, 
but the parametric decay heats the central part of the disc, 
and the corresponding pressure increase pushes matter away 
from the midplane (Fig. |l5| ). 

The same parametric decay may be found in Run 4 (Fig. 
[l^ ), which is done on the same grid as Run 3, but starting 



from a laminar state with only the epicyclic shear flow. In 
Run 3 the parametric decay occurs between orbits 58 and 
60 (Fig. but it is followed by an exponential damping 
caused by the turbulence itself (Fig. |^. The damping rate 
for Run 3 as estimated in Sect. 3.2 is based on a time interval 
after the end of the episode of parametric decay, and there- 
fore represents the turbulent damping rate. Also, in Runs 3 
and 4 we find an enhanced heating of the central regions of 
the disc and a resulting density reduction (Fig. [l^ ). 



4.2 Application to warped accretion discs 

We now investigate the implications of our results for the 
large-scale dynamics of a warped accretion disc. In linear 
theory we may estimate the amplitude of the epicyclic mo- 
tion as (cf. Eq. |l^) 

^ ~ A, (30) 

Cs Q!v 

where A is the dimensionless amplitude of the warp. This es- 
timate is valid for a thin and sufficiently viscou s disc, that is 
for H/r Qiv ^ 1 (Papaloizou & Pringle 1983). For observ- 



able warps in which A exceeds the aspect ratio of the disc, 
the epicyclic velocities are comparable to, or greater than, 
the sound speed, as is the case in our numerical simulations. 
Based on the high velocities one might expect shocks to ap- 
pear in the simulations, and shocks do appear in some global 
simulations that allow horizontal gradients in the epicyclic 
velocity (cf. Nelson & Papaloizou 1999). However, because 
of the local nature of our model we do not find any such 
gradients or shocks. 

The large-scale dynamics of a warped disc has been for- 
mulated by Pringle (1992) in terms of two effective kinematic 
viscosity coefficients: describes the radial transport of the 
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because the parallel component of angular momentum is 
transported mainly by these small-scale motions and mag- 
netic fields. However U2 oc instead of the intuitive result 
U2 oc av (Here we assume H/r av <C 1.) This is be- 
cause the perpendicular components of angular momentum 
are transported mainly by the systematic epicyclic motions 
induced by the warp, and these are proportional to as 
explained in Section 2.3. In effect, U2 is an effective viscosity 
coefficient removed by an additional level from the turbu- 
lent scales. If we generalize the linear theory of Papaloizou 
& Pringle (1983) to allow for an anisotropic small-scale ef- 
fective viscosity, we obtain v^jvx — l/(2QhQ!v)- 

The condition for a warp to appear in the accretion disc 
is set by the balance between the torque that is exciting the 
warp and the viscous torque, described by U2, that is flat- 
tening the disc. The warp-exciting torque may for instance 
be a radiation torque from the central radiation source. As- 
suming that accretion is responsible for all the radiation, 
the radiation torque will depend on the viscosity ui. The 
criterion for the warp to appear will t hen d epend on the ra- 
tio of viscosities 77 = V2/i^i- Pringle (1996) showed that an 



irradiation-driven warp will appear at radii 



Figure 16. {u'^)v (top) and («^)v (bottom) for Run 3 
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Figure 17. The density (p) (top) and internal energy (e) (bot- 
tom) as functions of 2 for Run 3 at t = 55.8 (solid line) and 
t = 60.9 (dashed line) orbital periods 



component of the angular momentum vector parallel to the 
tilt vector, while V2 describes the transport of the perpen- 
dicular components. The relation between and U2 and a is 
non-trivial but has been explained by Papaloizou & Pringle 
(1983) and Ogilvie (1999). Recall that ah and are ef- 
fective viscosity parameters that represent the transport of 
momentum by turbulent motions (and magnetic fields) on 
scales small compared to H. Now i/i oc Qh as expected. 



where i?sch is the Schwarzschild radius, and e — L/Mc? 
is the efficiency of the accretion process. We have shown 
that Oh ~ Qfv ^ 1- However, we emphasize again that this 
does not imply that « 1; on the contrary, we estimate 
that rj ~ l/(2Qhav) 3> 1. The high value for r) will make it 
difficult for a warp to appear unless the radiation torque can 
be amplified by an additional physical mechanism. One way 
to produce a stronger torque is if the irradiation is driving 
an outfiow from the disc (cf. Schandl & Meyer 1994). 

A similar damping mechanism may affect waves excited 
by Lense-Thirring precession in the inner part of the accre- 
tion disc around a spinning black hole. Numerical calcula- 
tions bi;_Markovic & Lamb (1998) and Armitage & Natara- 



jan (199E) show that these waves are damped rapidly unless 
1^2 <^ vi, which we find is not the case. (However, we note 
that the resonant enhancement of V2 will be reduced near 
the innermost stable circular orbit, because the epicyclic fre- 
quency deviates substantially from the orbital frequency). 
Likewise a high value of 1^2 will lead to a rapid alignment 
of the angular momentum vectors of a black hole and its 
surrounding accretion disc (cf. Natarajan & Pringle 1998). 

Some caution should be exercised when interpreting the 
results of this paper. Although the shearing box simulations 
in general have been successful in demonstrating the ap- 
pearance of turbulence with the right properties for driving 
accretion, they are in general producing uncomfortably low 
values of ah to describe for instance outbursting dwarf nova 
discs (e.g. Cannizzo, Wheeler & Polidan 1986). An under- 
estimate of ah and would lead to an overestimate of 77. 
In addition the parametric instability leads to an enhanced 
damping of the epicyclic motion. If this effect had been sus- 
tained it would have resulted in a smaller value for r). In a 
warped accretion disc the epicyclic motion is driven by the 
pressure gradients, which may maintain the velocities at a 
sufficient level for the parametric instability to operate con- 
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that the timescale for damping a warp in the accretion disc 
is much shorter than the usual viscous timescale. That the 
two alphas are within a factor of a few of each other is sur- 
prising, since the damping of the epicyclic motion may be 
attributed to the Reynolds stresses, while the accretion is 
mostly driven by the Maxwell stress. 

However, not all of the damping can be described as 
a simple viscous damping. At high amplitudes the epicyclic 
motion may also decay parametrically to inertial waves. The 
damping is much more efficient in the presence of this mech- 
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Figure 18. The upper plot shows {B"^), the average of the mag- 
netic energy, as a function of z for Runs 3 at 66.3To (sohd line) 
and at 66.IT0 (dashed hne). The lower plot shows —{BxBy) 
plotted the same way 
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tinuously. To address this question we intend to study forced 
oscillations in a future paper. 



4.3 The vertical structure of the accretion disc 

Our previous simulations of turbulence in a Keplerian shear- 
ing box have shown that the turbu lent a:t/-stresses are 



approximately constant with height (Brandenburg et al 



199(: ) rather than proportional to the p ressure as may have 
been expected from the a-prescription ( Shakura & Sunyaev 



197; ). We modified the vertical boundary conditions for this 
paper and added an epicyclic motion. The {BxBy)-stiess is 
still approximately independent of z or even increasing with 
1^1 for 1 2: 1 < H though, while at larger \z\ we may see the 
effects of the boundary conditions (Fig. II8I). The effect of 
the epicyclic motion is seemingly to limit the {BxBy)-stress 
in the surface layers to its value in the interior of the disc. 
The fact that the stresses decrease more slowly with z than 
the density results in a strong heating of the surface layers. 



5 CONCLUSIONS 

In this paper we have studied how the turbulence in an ac- 
cretion disc will damp an epicyclic motion, whose amplitude 
depends on the vertical coordinate z in the accretion disc. 
Such a motion could be set up b y a warp in the accretion 
disc (Papaloizou & Pringle 1983). We find that the typical 



damping timescale of the epicyclic motion is about 25 or- 
bital periods, which corresponds to = 0.006. This value 
is comparable to the traditional estimate of Oh that one gets 
from comparing the {pUxUy — BxBy}-stTess with the pres- 
sure. Both alphas are of the order of 0.01, which implies 
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